Simulation of Cu-Mg metallic glass: Thermodynamics and Structure 
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We have obtained effective medium theory (EMT) interatomic potential parameters suitable for 
studying Cu-Mg metallic glasses. We present thermodynamic and structural results from simulations 
of such glasses over a range of compositions. We have produced low-temperature configurations by 
cooling from the melt at as slow a rate as practical, using constant temperature and pressure 
molecular dynamics. During the cooling process we have carried out thermodynamic analyses based 
on the temperature dependence of the enthalpy and its derivative, the specific heat, from which the 
glass transition temperature may be determined. We have also carried out structural analyses using 
the radial distribution function (RDF) and common neighbor analysis (CNA) . Our analysis suggests 
that the splitting of the second peak, commonly associated with metallic glasses, in fact has little to 
do with the glass transition itself, but is simply a consequence of the narrowing of peaks associated 
with structural features present in the liquid state. In fact the splitting temperature for the Cu-Cu 
RDF is well above Tg. The CNA also highlights a strong similarity between the structure of the 
intermetallic alloys and the amorphous alloys of similar composition. We have also investigated the 
diffusivity in the supercooled regime. Its temperature dependence indicates fragile-liquid behavior, 
typical of binary metallic glasses. On the other hand, the relatively low specific heat jump of around 
1.5fcs/at. indicates apparent strong-liquid behavior, but this can be explained by the width of the 
transition due to the high cooling rates. 
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I. INTRODUCTION 

Metallic glassesii^ have generated considerable scien- 
tific interest since they were discovered 40 years ago, due 
to their unusual magnetic and mechanical properties, as 
well as wear and corrosion resistance, and their glass- 
forming ability per se. This interest has substantially in- 
creased since the discovery of the so-called bulk metallic 
glasses (BMGs) or bulk amorphous alloys, by Inoua^ and 
Johnson. i The ability to create samples with thicknesses 
in the mm or cm range, rather than /im thick ribbons, 
greatly increases the applicability of the materials, as well 
as the range of measurements that can be performed on 
them. This is particularly true in the case of mechani- 
cal testing, and recently measurements of properties such 
as fracture toughness, fracture morphology and crack-tip 
plasticity have been madei^*SiL& 

The mechanisms of plastic deformation are of partic- 
ular interest in metallic glasses in view of the fact that 
there are no obvious topological defects which might play 
a role analogous to crystal-dislocations, allowing slip to 
take place in small increments. Thus metallic glasses 
tend to have very high flow stresses^ A complete un- 
derstanding of plastic deformation must include the fol- 
lowing two parts: (i) detailed knowledge of the elemen- 
tary events that constitute plastic flow and (ii) a practi- 
cal continuum theory which uses this knowledge to make 



predictions of macroscopic behavior (a recent such the- 
ory is the so-called Shear Transformation Zone (STZ) 
theorySJ^i) . The motivation for the present work is a 
desire to tackle item (i) using the tools of modern ma- 
terials simulations, specifically: realistic potentials, sys- 
tem sizes as large as feasible and necessary, and sophis- 
ticated analysis and visualization techniques. The first 
step, addressed in this paper, is to create appropriate in- 
teratomic potentials, generate glassy configurations, and 
study the thermodynamics and structure of the system, 
in order to understand it as a glass-forming one. Sim- 
ulations of mechanical properties will be presented in 
future publications. The phrase "realistic potentials" 
refers to contemporary potentials commonly used for 
metals, including effective medium or embedded atom- 
type potentials, or pseudopotential-based pair potentials, 
as opposed to Lennard-Jones potentials, which are com- 
monly used (with two components) to model metallic 
glassesi -^^i^^i^^i^'^i^^i^^ Such potentials are especially use- 
ful because they allow quantitative comparison with ex- 
periments of properties such as glass transition temper- 
ature, and, later, mechanical properties. 

In this paper we present molecular dynamics simula- 
tions of the binary alloy Cua:Mgi_a;. Mg-based BMGs 
such as Mg6oCu3oYio2ii2ii2iiS are of interest because 
their weight is low, being dominated by Mg, but their 
strength can be comparable to high-strength steel. We 
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have chosen to study the binary alloy because (i) it is 
simpler to optimize a potential for two species than for 
three and (ii) it is easier to study dependence on a sin- 
gle composition-parameter than on two. Our intent to 
use realistic potentials necessitates an attempt to create 
as realistic a glass as possible with those potentials. It 
is thus important to characterize the system as a glass- 
forming and alloying one as completely as possible. 

The Cu-Mg equilibrium phase diagram is shown in 
Fig. ^ Experimentally it forms a glass over a range of 
compositions from 9-42 at. % (complete glass formation 
over 12-22 %, which includes the eutectic composition 
14.5 It is not a BMG, since it can only be formed 

by melt-spinning at high cooling rates. The cooling rates 
in the simulations are necessarily even higher and allow 
glassy configurations to be created over almost the entire 
range of compositions. It is worth studying the experi- 
mentally inaccessible states as part of the process of de- 
tecting trends in material properties as a function of com- 
position; it is the crystal-nucleation timescale, lying be- 
tween the simulation and experimental timescales, which 
makes the difference between crystalline and amorphous 
phases — if just a few orders of magnitude gain in cool- 
ing rate could be experimentally realized, there is reason 
to believe these states would be as stable as the actual 
glassy configurations currently realizable by experiment. 

Because the crystallization rates are high, there are 
limited experimental measurements of the thermody- 
namic properties of Mg-Cu glasses, and thus it is of inter- 
est to study these in the simulations before moving on to 
mechanical properties. In the process we find some inter- 
esting results regarding structural changes in the super 
cooled regime (steady growth of icosahedral order and 
evidence of restructuring thermodynamics). Addition- 
ally we make some observations on the question of the 
fragility of this system. The next section will discuss 
some aspects of the theory of glass formation in alloys, 
as applied to the Cu-Mg system. Section ITTll will discuss 
simulation methods, including the fitting of the inter- 
atomic potential. Sections IIVI and Ivl discuss characteri- 
zation of the glass transition and of structural properties 
respectively. The last section is the discussion. 



II. GLASS FORMATION IN THE Mg-Cu 
SYSTEM 

One approach to the theory of metallic glass for- 
mation is based on pseudopotential-derived interatomic 
pair potentials, ^^'^^ and emphasizes the coincidence of 
bond lengths with potential minima. We will not be 
using such potentials; in fact many aspects of glass 
formation are purely geometrical (packing of spheres) 
and phase-energetioSi (comparison with competing crys- 
talline phases). Frank and KaspeoSSi^ pointed out that 
many complex intermetallic structures can be understood 
in terms of tetrahedral close-packing of spheres. Exam- 
ples of so-called Frank-Kasper (FK) phases include the 



Laves phases (C14, 015, 036) and /i, x f^nd cr phases. 
The high packing fraction and coordination numbers sug- 
gest directional bonding does not play a role. The closest 
packing of spheres of equal size is achieved with a tetra- 
hedron (79%), but tetrahedra cannot fill space — the best 
one can do is make an icosahedron out of twenty slightly 
distorted tetrahedra, but this cannot be repeated period- 
ically, so in crystals one has the fee (e.g. Ou, a = 3.6lA) 
and hep (e.g. Mg, a = 3.2lA, c = 5.2lA) structures, 
with 74% packing. 

In the Ou-Mg system there is indeed a Laves phase, 
Ou2Mg. This is not surprising given that the ideal Laves 
packing is achieved with a radius ratio of 1.225 fRef. l2^ 
p. 59), which is close to that of Mg and Ou (1.256 using 
the Goldschmidt radii, based on nearest neighbor dis- 
tances of the pure metals). This phase is quite stable 
simply because having a majority of smaller atoms allows 
a greater packing fraction. On the other hand Mg20u, 
with the larger atoms in the majority, is not as stable 
an alloyi2i Mg-Ou is in a class of metallic glass formers 
which include simple metal-transition metal binary alloys 
and are characterized by a Laves phase when the small 
atom (Ou) is in the majority, and a glass when the larger 
atom is. In Ou2Mg the Ou atoms have 0N12 icosahe- 
dral coordination and the Mg atoms are 16-coordinated, 
surrounded by so-called Frank-Kasper 16-hedra (more 
specifically, Friauf polyhedra). 

Glass formation in a binary alloy appears to be fa- 
vored by the same criteria that favor the formation of 
FK phases: large, negative heats of formation, non- 
directionality of bonding, and a tendency to maximize 
packing fraction. In general one finds that for composi- 
tions between intermetallics (for example near eutectics) , 
where the equilibrium phase diagram shows a two-phase 
mixture, the amorphous phase is more stable than any 
single crystalline phase. In the region of the phase dia- 
gram where FK phases appear, glass formation typically 
loses out in the competition experimentally, presumably 
because the nucleation of the Laves phase is rather easy. 
In the Ou-Mg system, the region of experimental glass 
formation is on the Mg-rich side, where the competing 
crystalline phase, Mg20u, is quite complex (48 atoms in 
the unit cell). 



III. SIMULATION METHODS 

A. Potentials 

The interatomic potential we use is the effective 
medium theory (EMT),^^'^^ fit to data obtained from 
density functional theory (DFT) calculations and exper- 
iment. This has previously been applied to fee metals, 
in particular late transition and noble metals and has 
been of great use in studying mechanical properties of 
crystalline metalsi^SiSS As Mg crystallizes in hep with an 
almost ideal c/a-ratio of 1.624 (ideal is ^8/3 = 1.633), 
indicating little directional bonding, we might expect it 
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FIG. 1: Equilibrium phase diagram for Mg-Cu (from Ref . l25h 
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TABLE I: EMT parameters for Cu and Mg, in units derived 
from eV and A. 



to be reasonably well described by an appropriately op- 
timized EMT potential. 

EMT uses seven parameters for each element. A set of 
parameters for Cu exists but these have been optimized 
for simulations of pure, crystalline Cu (where for exam- 
ple particular attention was paid to the stacking fault 
energy, which is of no concern in amorphous materials). 
For the amorphous alloys, it is important that the forma- 
tion energies are reasonable, in particular that they are 
negative (otherwise the system will simply separate into 
regions of pure Cu and regions of pure Mg) . 

Thus we have (re-)fit the parameters of both elements, 
taking into account basic properties such as lattice con- 
stants, cohesive energies and elastic constants of the pure 
elements, as well as the formation energies of the two in- 
termetallic compounds, Mg2Cu and Cu2Mg. Due to the 
near-ideal hep packing of Mg, its structure differs from 
fee only at the second neighbor level. For simplicity, and 
because the EMT potential is formulated in terms of fee 
packing, we used calculated properties of fee Mg in the 
fitting, except that the cohesive energy was corrected us- 
ing the experimental hep value and the calculated fcc- 
hcp difference (23mey/atom), calculated differences in 
cohesive energy being expected to be more accurate than 
calculated cohesive energies themselves. 
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Cu2Mg-Aii" 


-0.159 
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TABLE II: Properties used in the fitting: the values speci- 
fied (from DFT/experiment) and the values according to the 
optimized potential. B = bulk modulus, a — lattice constant. 

The optimized EMT coefficients are shown in table |2 
and the target and fitted values of the fitting proper- 
ties are shown in table |n| Note that for the orthorhom- 
bic Mg2Cu, the experimental b/a and c/a were used, as 
well as the experimental values of the internal coordi- 
nates. The alloy formation energies are well-represented. 
Unlike pair potentials based upon pseudopotentials, the 
present form of the EMT potential^ does not incorpo- 
rate the Friedel oscillations, and the idea that stability of 
intermetallic compounds is determined by the matching 
of minima of pair potentials to interatomic distances^ 
does not play a role; the fact that EMT parameters can 
be chosen to give the correct formation energies of the 
intermetallic compounds appears to be most important. 



B. Molecular dynamics 

We simulated the cooling of systems of 2048 atoms 
from the liquid state (above the melting point) down to 
zero temperature. The compositions ranged from pure 
Mg to pure Cu, and are labeled by the percentage of Cu. 
For most simulations we used 21 compositions, increas- 
ing in steps of 5% from (pure Mg). The initial state 
was an fee lattice with the sites occupied at random by 
Cu or Mg atoms in accordance with the overall composi- 
tion. There was no initial heating phase; the first stage 
in the cooling run set the temperature to a value well 
above the melting point (values ranged from 1392 K for 
Mg {Tm = 923K) to 1857 K for Cu (T™ = 1358K)), 
making the crystal melt immediately. Two rates of cool- 
ing were used; differing in the amount of simulation time 
at each temperature. Cooling took place in steps of 
35K; the procedure at each temperature stage was as 
follows: (i) a small number of steps, corresponding to 
0.6ps (the MD time-step was 2/s), of constant-volume 
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Langevin thermalization was carried out in order to ap- 
proximately thermalize the system to the new tempera- 
ture; (ii) the dynamics was switched to constant-pressure 
(NPT) dynamics and the system was simulated for an 
initial equilibration time of 6ps/12ps; (iii) the system 
was simulated for a longer time A0ps/120ps during which 
thermal averages of various quantities of interest were 
taken. This time also contributed to the equilibration of 
the system. The overall cooling rates were thus close to 
0.72K/ps (7.2 X lO^^K/s) and 0.25K/ps (2.5 x lO^if/s). 
The NPT dynamics used was a combination of Nose- 
Hoover and ParrincUo-Rahman dynamics, proposed by 
Melchionnai^'^^''^^''''^ We turned off shearing, allowing 
only volume fluctuations, because the liquid state can- 
not support a shear stress and fluctuations in the peri- 
odic box sometimes lead to extreme angles between box 
vectors and thus problems with the neighbor-locating al- 
gorithm. The pressure was zero or a small positive value 
(this was necessary in some cases when the initial tem- 
perature was above the boiling point of pure Mg). For 
each cooling rate the simulations were run twice with dif- 
ferent random number seeds (affecting the distribution of 
species in the initial lattice and the Langevin dynamics 
used when the temperature is changed; the NPT dynam- 
ics does not use random numbers). 

During the averaging period, the pressure, volume, ki- 
netic and potential energies were recorded and averaged. 
For the purposes of structural analyses so too was the 
radial distribution function (RDF), both total and sepa- 
rated into contributions from Mg-Mg, Cu-Cu and Mg-Cu. 
At the end of the averaging time the current configura- 
tion was saved, as well as a configuration obtained from 
it by direct minimization (quenching) using the MDmin 
minimization algorithm. At a later time the saved con- 
figurations from selected temperatures were used for fur- 
ther simulation at that temperature to gather further dy- 
namical and structural information such as diffusion con- 
stants and thermally averaged common neighbor analysis 
(CNA). 

Our cooling rates are as slow as in other recent sim- 
ulations of amorphous metalsj^^iSSiSSiSl but they are of 
course larger than experimental rates by several orders 
of magnitude. In order to check that our results are not 
significantly affected by this difference, we have cooled 
one composition, 15% Cu at several faster rates and one 
slower one. Fig. [21 shows Tg and the enthalpy at T = 
for these runs. The methods of calculating Tg are ex- 
plained in the next section; only one (intercept) could be 
used for the very fast runs. It is pretty clear that for 
the cooling rates used in the main simulation, the depen- 
dence of Tg on cooling rate has become smaller than the 
uncertainty in determining Tg. The enthalpy shows a def- 
inite slope still at the lowest cooling rate, amounting to 
about ImeV per order of magnitude cooling rate, which 
is rather small; also one would expect the curve to flatten 
out more at even smaller rates. The one significant differ- 
ence we notice is that crystallization at the Cu-rich end 
happens at lower Cu-concentrations for slower cooling: 
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FIG. 2: Cooling rate dependence of Tg for 15% Cu system. 
Solid symbols, Tg from the maximum rate of change of Cp; 
open symbols, Tg by intercept method. Arrows indicate the 
cooling rates used for the main simulations. Inset: enthalpy 
of system at end of cooling run (T = 0) . 
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FIG. 3: Specific heat versus temperature for 15% Cu system, 
cooling rate 0.72K/ps. Dashed lines, values from two separate 
cooling runs, displaced ±0.2 for clarity. Solid line, average of 
these two. Solid vertical line, Tg from maximum slope of 
specific heat; solid dashed line, Tg from intercept method. 
Inset: enthalpy (average of the two runs) versus temperature. 
Dotted lines show the extrapolated straight-line fits from the 
intercept method. 



the 90% Cu system crystallizes in one run at 0.25K/ps 
but not at all at 0.72K/ps. 



IV. GLASS TRANSITION 

We see glass transitions in almost all compositions, the 
exceptions being the pure elements and 95% Cu, which 
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crystallize in fcc/hcp structures (also 90% Cu in one out 
of two runs at Q.2^K /ps). The first evidence that a glass 
transition takes place upon cooling appears in the en- 
thalpy versus temperature curve, which shows a change 
in slope (inset in Figl^J. This suggests a way to deter- 
mine Tg by breaking the curve into two pieces, fitting a 
straight line to each, and intersecting the two lines ob- 
tained. We call this the "intercept-method". It turns 
out that this tends to underestimate Tg as can be seen 
by looking at the derivative of the enthalpy, the specific 
heat (Fig 13), obtained from centered differences. The 
Tg ends up at the leftmost part of the steep part of the 
curve, whereas one would expect any reasonable defini- 
tion of Tg to be roughly in the center of the transition 
region (defined as the steep part). Thus we compute Tg 
as the temperature at which the specific heat is chang- 
ing fastest by taking derivatives again and simply choos- 
ing the maximum. This method necessarily yields a Tg 
equal to one of the simulation temperatures, but since 
the transition region is a few times wider (150 — 200 A') 
than the temperature step in the simulation, one cannot 
expect to do better (experimentally one sees widths of 
some tens of K, see for example Ref . Is^ . In cases where 
we have two different enthalpy curves for the same system 
cooled identically but from different starting configura- 
tions we average the two enthalpy curves before applying 
the method, as this gives a smoother cp curve. 

The Tg we get for 15% Cu is ibQK which is remark- 
ably similar to the experimental value of iSQK reported 
by Sommer et al^ In runs where crystallization took 
place, a large spike in the specific heat appeared, corre- 
sponding to a step or latent heat in the enthalpy curve. 
Before looking at the composition dependence of Tg, we 
notice that the temperature dependence of cp is quite 
similar in form to experimental specific heat curves of 
Zr4i.2Tii3.8Cui2.5Niio.oBe22.5 reported by Busch et alu^^ 
and of fluorozirconate and tellurite glasses reported by 
Lin and Navrotsky)^2il£ there is an increase in specific 
heat in the supercooled liquid region compared to the 
high-temperature liquid region. For the tellurites, Lin 
and Navrotsky identified the source of this as specific 
structural rearrangements that take place in the liquid 
prior to the glass transition. We will see in the next sec- 
tion what evidence there is for structural rearrangements 
in the Cu-Mg supercooled liquid. 

Fig.^shows Tg and Acp, the heat capacity jump (ob- 
tained by roughly determining the transition region as 
the peak in the derivative of Cp and taking the difference 
of Cp on either side of the peak) for different composi- 
tions and cooling rates. Tg rises roughly linearly with 
increasing fraction of Cu, which presumably reflects a 
general increase in energy scale as we go from the weakly 
cohesive (low melting point) Mg to the more strongly 
cohesive Cu. The fluctuations towards the Cu-rich end 
are due to the mid-point method's difficulty handling the 
somewhat less clean cp-data there. The fluctuations in 
Acp are also due to the imperfect cp-data. Nevertheless, 
it seems clear that Acp has the value of roughly X.bks 
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FIG. 4: Tg and Tput (upper panel) and Acp (lower panel). 
Squares: Tg, Acp at 0.72K/ps; diamonds at 0.25K/ps. Tri- 
angles: rspiit(Cu) at 0.72K/ps. 

Enthalpy at K 



-10 



-15 



o o 
o oo 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

fraction Cu 



FIG. 5: Diamonds: formation enthalpy per atom in final zero- 
temperature glassy state. Dotted line: formation enthalpy per 
atom of corresponding (in general two-phase) crystal. 



per atom, independent of concentration. This is a rela- 
tively small amount, which is typical of so-called "strong" 
glass-formers, which include most BMGs4i In particular 
the Mg65Cu25Yio shows a jump of the same order (ac- 
tually 2fcB/at.).^^ However, we should be careful about 
inferring strong-liquid behavior from this measurement; 
binary alloys typically are not strong glass-formers^^ and 
below we shall see evidence of fragile-liquid behavior in 
the diffusivity. The apparent small jump of Cp may be 
a consequence of the width of the transition. 

As a partial means of determining how good, mean- 
ing how stable or well-annealed, the final configurations 
are, we consider their enthalpies. We have seen already 
how the final enthalpy depends on cooling rate (Fig. |2J); 
we now compare to the equilibrium phases, for different 
compositions. Fig.|Slshows the formation enthalpies as a 
function of composition. The formation enthalpy is the 
enthalpy minus the appropriate linear combination of the 
pure elements' enthalpies. The appropriate quantity to 
compare to, also shown in Fig. [S] is the formation en- 
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configurations for further simulations in which diffusion 
constants for the two atomic species were measured. An 
Arhcnnius plot for the 15% Cu composition is shown in 
Fig There is a clear indication of a transition near 
1000/T - 3K-^, corresponding to T ~ 330i^, which is 
consistent with the Tg — 350K obtained from the specific 
heat. For each composition for which diffusion constants 
were measured, we have fitted the high temperature part 
of the data to the Vogel-Fulcher (VF) law 
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FIG. 6: Diffusion constants in 15% Cu. Squares: Mg; dia- -C = ^oexp( ^ _ ^ ) (1) 

monds: Cu. Dotted line: VF fit to high temperature Mg data. ^ 
Dashed line: Arrhenius fit to low temperature Mg data. 
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FIG. 7: Upper panel: fragility parameter D* of Vogel-Fulcher 
fits to Cu diffusion constants at selected compositions. Lower 
panel: location of apparent divergence To from the same fits. 



thalpy of the corresponding crystalline phase, which in 
general is a two-phase mixture (so, e.g., between 33% 
Cu and 66% Cu it is an appropriate weighting of the 
formation enthalpies of Cu2Mg and Mg2Cu). We notice 
that the glass formation enthalpy follows quite closely 
the crystalline one, being 1-4 kJ/mol higher (the excep- 
tions being at 0, 95% and 100% Cu where the system 
did in fact crystallize). This is quite small and typical 
of easy glass formersiSii^ For the 15%Cu composition 
the value 4:.2kJ/mol was reported by Sommer et al. for 
the transformation enthalpy from the crystalline to the 
amorphous state, which is slightly higher than our value 
of 3.53k J/mol — that is, our glass at this composition ap- 
pears to be a little too stable compared with experiment. 
This kind of discrepancy can only be due to limitations 
of the interatomic potential, and not to the high cooling 
rate. This gives us further confidence that we have cre- 
ated glassy structures which are more or less as stable as 
they can be. 

For selected compositions and selected temperatures, 
configurations from the cooling runs were used as initial 



where Tq is the location of the apparent singularity and 
D* is the so-called fragility parameter. In Fig[3we show 
D* and Tg obtained from fits of the Cu diffusion con- 
stants to the VF law (the Mg values are very similar, the 
differences being very small compared to the differences 
from composition to composition) . There is a reasonably 
clear trend towards decreasing D* and increasing Tq as 
the fraction of Cu increases. High values D* are asso- 
ciated with strong glass-formers, the archetypal case of 
Si02 having D* = 100. Bulk metallic glasses are con- 
sidered strong^^ with D* ^ 20. So-called fragile glasses 
have D* around 2. From our diffusion data we get low 
fragility parameters, in the range 2-4, indicating that the 
Mg-Cu glasses are somewhat fragile. This is consistent 
with the experimental fact that this is not in fact a bulk 
metallic glass. The Tq values increase as the fragility de- 
creases, so that the apparent singularity approaches the 
actual glass transition temperature. These trends, re- 
flecting greater fragility (decreasing D*) with increasing 
Cu composition, are also consistent with the fact that 
experimentally, amorphous Mg-Cu can only be made at 
all for Mg-rich compositions, since strong liquids tend to 
be robust against crystallization (in a strong glass for- 
mer the melt viscosity is high, making the kinetics slow) . 
Thus, our diffusion results put the binary alloy Mg-Cu 
at the fragile end. This seems to contradict the sugges- 
tion of strong-liquid behavior from the specific heat data. 
The small ACp, may have a simple explanation, however, 
namely that it has been reduced due to the broadening of 
the transition in the simulations compared to what one 
would expect experimentally. This broadening implies 
that a certain amount of restructuring, that at slower 
cooling rates would take place above Tg, in the simula- 
tion takes place during and below Tg. The net enthalpy 
change (the area under the Cp curve) is more or less the 
same, so the height of the curve above Tg must be reduced 
to compensate. Thus we can assert that the simulations 
are consistent with Mg-Cu being a fragile glass-former, 
like most binary alloys. 
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FIG. 8: RDFs, Mg-Mg and Cu-Cu. Left panels: 15% Cu; 
right panels: 85% Cu; upper panels: T = Tout; lower panels: 
T — Tg\ inset on bottom left panel: combined RDF (Mg-Mg, 
Cu-Cu, Mg-Cu). 



V. STRUCTURAL ANALYSIS 
A. Radial distribution function 

Fig.lHlshows the partial RDFs 5Mg-Mg(?') and gcu-Cu{i') 
for two compositions at two temperatures. At the higher 
temperature, which is the eutectic temperature for the 
corresponding region of the phase diagram, the system 
is expected to be in equilibrium, and the RDFs have the 
normal structure of a liquid, with nearest neighbor dis- 
tances of 3.lA for Mg and 2.6A for Cu, which are close 
to their values in the bulk crystal phases of the pure ele- 
ments. The lower panels in Fig. [S] show the RDFs at the 
respective Tg for each composition. At 15% Cu, the first 
peak is prominent for both RDFs. In the Cu-rich alloy 
on the other hand, the first Mg-Mg peak is significantly 
suppressed, indicating the Mg atoms are not particularly 
likely to be found next to each other. This is not surpris- 
ing since we expect Mg-Mg bonds to be weak compared 
to both Cu-Cu and Mg-Cu bonds, given the cohesive en- 
ergies of the pure elements and the intermetallic com- 
pounds. 

We can see a distinct splitting of the second peak in 
ffCu-Cu in both compositions. The splitting occurs also 
for 5Mg-Mg, but at lower temperatures (here it is also ob- 
scured, particularly in the Cu-rich compositions, by the 
fact that the first sub-peak is significantly higher than the 
second, which thus appears as a shoulder on the high-side 
of the first). Such a splitting is commonly associated with 
the glass transition, but we can see here that the splitting 
is already well developed at Tg for gcu-Cu and in fact it 
first occurs well above Tg. Fig. ^ shows the temperature 
^spiit at which this occurs, determined in a somewhat ar- 
bitrary manner by visual inspection of the RDFs for dif- 
ferent temperatures, as a function of composition. The 
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FIG. 9: Partial and total coordination numbers as indicated. 
ZaB means the average number of neighboring B atoms that 
an A atom has. Crosses are Spaepen-Cargill short-range 
chemical order parameter rjAB determined from the coordi- 
nation numbers. 



dependence on composition is rather less than that of Tg, 
and in fact it appears that the splitting is not related to 
the glass transition in a direct way. Note that what is 
typically observed experimentally — the combined RDF, 
which averages over the different components — does not 
show the splitting, because the location of the second 
peak differs for different components and the effect is 
washed out (see inset Fig. |S1 in fact, for Mg-rich compo- 
sitions, the total RDF has a split first peak due to the 
difference in location between Cu-Cu and Mg-Mg first 
peaks). 



B. Coordination numbers 

By integrating the RDFs appropriately^ we can de- 
termine the partial, total and average coordination num- 
bers, Zab, Za and Z. These are shown in Fig.|51 for the 
zero temperature RDFs from the runs with the higher 
cooling rate (0.72/^/^5). We have checked that virtu- 
ally identical results are obtained with the lower cooling 
rate. The average coordination number is quite indepen- 
dent of composition, Z = 12.91 ±0.17. The coordination 
number of Mg, Zjv/g, is always higher than the total Z^ 
and Zcu always lower. Both rise as the fraction of Cu 
increases (their average does not because it is weighted 
by the concentrations). From the coordination numbers 
we can calculate the Spaepen-Cargill short-range order 
parameter- 



,43 



= ZAB/Z*Ag — 1 
^*AB = CbZaZb / Z 



(2) 
(3) 



Here cb is the concentration of B atoms (which we take 
as Cu) . A positive value of t^ab indicates a tendency for 
more unlike bonds than would be expected in an alloy 
which is completely chemically disordered. Fig. |5| shows 
rjAB ; it is definitely positive throughout the glassy range 
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FIG. 10: Common Neighbor Analysis (CNA) of first peak of 
partial RDFs for 50%Cu glass. Bottom right panel; RDF 
(solid line) contributions from 555, 544, 433, 421 and 422 
pairs (dotted lines), and the sum of these (dashed line). Other 
panels: number of neighbors of specified type (e.g. of a Cu 
atom which are Cu and make a 555 pair) as a function of 
temperature. Squares: 555, diamonds: 544, triangles: 433, 
plus: 421; cross: 422. 

of compositions, with an apparent maximum near the 
middle of the range. However while the maximum value 
r]AB can ever take is unity, the maximum for a particu- 
lar composition is somewhat less, and values should be 
normalized by the maximum before comparing different 
compositions. We have not done this since it is not clear 
what tf^g^ is in an amorphous system, particularly in 
the regime of a;^ — xb (see Ref . 143) . 

C. Common Neighbor Analysis 

To obtain more detailed information about the lo- 
cal atomic structure we use Common Neighbor Analysis 
(CNA) 4^*^ The analysis assigns three indices to every 
pair of atoms, thus allowing a decomposition of the RDF 
into contributions from different types distinguished by 
their CNA indices. The first index is the number of neigh- 
bors that the two given atoms have in common; the sec- 
ond is the number of bonds among those neighbors, and 
the third is the size of the largest bonded cluster within 
the common neighbors (this last differs from the original 
definition,-^ but agrees in all the cases of interest and 
is less ambiguous). The cutoff for two atoms to be con- 
sidered "neighbors" or "bonded" is the position of the 
first minimum in the appropriate RDF. Note that the 
separation of the two atoms to whom the indices are as- 
signed can be anything up to twice the nearest neighbor 
distance — beyond this, they cannot have any neighbors 
in common. Several groups have presented CNA analyses 
of the structure of metallic glasses i^^i^SilS-^SiL^ These 
all reported similar results: the first peak of the RDF is 



composed mostly of 555, 544 and 433 pairs, and the sec- 
ond peak is composed mostly of 333, 211 and 100 pairs. 
555 pairs are associated with icosahedral order: in a per- 
fect icosahedron the central atom makes a 555 pair with 
each of its 12 neighbors. 544 and 433 pairs are formed 
when one or more bonds between the outer atoms of an 
icosahedron are broken. 333, 211 and 100 pairs can also 
be associated with various pairs within a perfect icosahe- 
dron. Furthermore, the 333 and 211 pairs of the second 
peak combine to form the first sub-peak and the 100 pairs 
make the second sub-peak, when the second peak splits. 

In these papers the CNA was always performed on 
quenched configurations, obtained by rapid minimiza- 
tion to local minima from finite temperature configu- 
rations; this is preferable to doing the analysis on an 
instantaneous configuration at finite temperature, since 
the distortions caused by thermal fluctuations would in 
that case obscure the "inherent structure" . Changes in 
the structure were correlated with the temperature from 
which the quench was made. In our analysis, we have 
taken an alternative approach to dealing with thermal 
fluctuations and have computed the full thermal averages 
of the contributions to the RDFs from pairs of different 
types. Analogous to the RDF which is itself a thermal av- 
erage, we thus obtain a "radial distribution function" for 
Cu-Cu/Mg-Mg/Mg-Cu pairs of type 555, 333, etc., which 
is in fact an exact decomposition of the full RDF for the 
given species-pair. These averages were computed dur- 
ing the same runs as the diffusion constants, with starting 
configurations taken from the cooling runs at 0.72K/ps. 
The CNA partial RDFs were computed every 10th ma- 
jor time step (starting with the 20th — after which it was 
assumed the full RDF had converged sufficiently to read 
the position of the first minima). The CNA partial RDFs 
each consist of a single peak from which quantities such as 
peak position, height and width can easily be extracted. 
Also computed is actual (average) number of pairs associ- 
ated with such a peak, obtained by integrating the RDF 
against 47rr^ times an appropriate density. Furthermore 
we can see directly how these CNA-RDFs sum to give 
the full RDF for a given species-pair. 

Fig. IIUI shows the numbers of pairs in the sub- 
components of the first peak of the RDFs. We see the 
same broad picture described above, in terms of the roles 
played by 555, 544, 333, etc. pairs. This should not be 
surprising since as we shall see later icosahedral order is a 
dominating feature of the intermetallic alloys. In partic- 
ular the number of 555 pairs grows more or less linearly 
as the temperature decreases from 1200 K to the glass 
transition temperature, beyond which it continues to in- 
crease, albeit with slightly smaller slope. The number of 
421 and 422 pairs, associated with crystalline hep and 
fee order, is very small at all temperatures. The bottom 
right panel of Fig. 1 101 shows the decomposition of the first 
peak of the RDF into contributions from the five listed 
pair-types. The difference between the solid line (full 
RDF) and the dashed line (sum of the five listed pair- 
types) indicates that other types make up a noticeable 
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FIG. 11: CNA of second peak of partial RDFs for 50%Cu 
glass. Bottom right panel: RDF (solid line) contributions 
from 333, 211 and 100 pairs (dotted lines), and the sum of 
these three (dashed line) . Other panels: number of neighbors 
of specified type (e.g. of a Cu atom which are Cu and make 
a 333 pair) as a function of temperature. Squares: 333, dia- 
monds: 211, triangles: 100. In the Mg-Cu case, the number 
refers to Cu neighbors of Mg atoms. 



fraction. This were found to include small amounts of 
311, 322, 666, 533, and 532 pairs. At high temperatures 
when the number of 555 pairs is low, all of these types 
of pairs, and some others not mentioned, contribute in 
small amounts to make up the full coordination num- 
bers. Thus the picture we have of the liquid structure at 
high temperatures is one of many (we have seen up to 15 
different CNA types for nearest neighbor pairs) different 
local structures constantly being created and destroyed, 
and all contributing a little bit to the thermal average. As 
the temperature cools, a pair of atoms is more and more 
likely to be found as a 555 pair. This is independent of 
what species the two atoms are, and of the compositions. 

Fig.llllshows a similar analysis of the second peak. We 
see what others have found previously, that it is mostly 
made up of the 333, 211 and 100, and the first two making 
up the first sub-peak and the latter the second sub-peak. 
In fact there is only a small difference between the sum 
of these three contributions and the full CNA, which ap- 
pears on the shorter-distance side of the peak. This small 
difference is due to 455, 444 and 322 pairs, which mainly 
occupy the region between first and second neighbor dis- 
tances. At the highest temperatures (not shown), these 
last three pairs make up a somewhat larger contribution, 
and are more clearly part of the second (main) peak, but 
the 333, 211 and 100 arc still dominant. The numbers 
of pairs associated with these CNA-types change rela- 
tively little with temperature: increases by about 
30% during cooling; iV2ii decreases by the approximately 
same amount, leaving their sum constant (A^ioo is in- 
volved in this only to a small extent). It seems that 
211 pairs are being transmuted to 333 pairs as the sys- 



tem cools. In the last section we saw that the specific 
heat of the supercooled liquid is higher than the high 
temperature liquid, and noted how such behavior in ex- 
periments, termed "restructuring thermodynamics" , has 
been associated with structural rearrangements that take 
place during cooling. In our system it is natural to as- 
sume that the rearrangements identified by CNA analysis 
in this and the preceding paragraph are responsible for 
the increased specific heat. 



D. Comparison with ordered structures 

At this point it is interesting to see what a common 
neighbor analysis of the intermetallic alloys Cu2Mg and 
Mg2Cu yields. The results are displayed in table IIIII 
along with the partial and total coordination numbers. 
The numbers of different CNA types could be separated 
further by the species of the second atom, but the table 
has already enough numbers. We see a distinct preva- 
lence of nearest-neighbor pairs of type 555 — almost all 
nearest neighbors are of this type, the rest being 444 and 
666. It is impossible to have every nearest-neighbor pair 
being of the 555 type in a crystal, but it certainly seems 
that the crystal structures here are trying to maximize 
the number of 555 neighbors. Now, "icosahedral order" 
strictly refers to having coordination number 12, all 555; 
however since in a binary alloy with a distinct size dif- 
ference this coordination number is only achieved for the 
smaller atom, and only in a certain composition range, 
strict icosahedral order cannot be attained, but we still 
choose to refer to a high number of 555 pairs as repre- 
senting "icosahedral order" . 

Of the second neighbor pairs only a few are of type 
211, most being 333 and 100. This is also consistent 
with icosahedral order: 333 pairs can be associated with 
pairs of tetrahedra which share a face, such as adjacent 
tetrahedra in a perfect icosahedron (or in the 555 struc- 
ture, in view of our generalized sense of "icosahedral"). 
211 pairs differ from 333 pairs by the removal of one 
of the common neighbors. It can be supposed that the 
211 pairs are defects of the icosahedral structure, just as 
544 and 433 pairs are, and thus that one would expect 
fewer of them relative to the number of 333 pairs in a 
more perfectly icosahedral structure. This is consistent 
with the fact that the numbers of 211 pairs decreases 
as temperature decreases in the glassy systems. In the 
same table are shown corresponding figures for the amor- 
phous alloys of closest composition to the intermetallics, 
except that the numbers of 555, 544 and 433 pairs have 
been combined under the 555 column. The numbers for 
the amorphous alloys agree with those from the corre- 
sponding crystalline phase to within 20 percent in most 
cases, the biggest difference being the reduced number of 
333 pairs, compensated more or less by the increase in 
211 pairs. If we were to combine the 333 and 211 fig- 
ures, like we have the 555, 544 and 433 ones, we would 
see that the structures in the amorphous and crystalline 
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TABLE III: CNA figures for nearest (Z) and next-nearest {N) neighbor pairs, for the intermetallic alloys Cu2Mg and Mg2Cu 
and amorphous alloys Mgo.35Cuo.65 and Mgo.50Cuo.50- Za and Zab are partial and total coordination numbers. For the 
amorphous alloys the figure under ^a555 represents a sum over 555, 544 and 433 pairs. 

on the short side. Thus the sphtting of the second peak 
does not itself indicate any kind of structural transition. 
It merely follows from the fact that the structure at this 
length scale (second neighbor distance) associated with 
the liquid state remains as one cools into the glass state, 
and the narrowing of peaks which is to be expected as 
thermal motion decreases. 



VI. DISCUSSION AND SUMMARY 

Our main intent in simulating the cooling of Mg-Cu al- 
loys has been to generate glassy configurations that can 
be considered realistic enough for simulations investigat- 
ing mechanical properties. Mg-Cu is a first step towards 
the more technologically interesting material Mg-Cu-Y. 
In order to assess the realism of the simulations we have 
studied various aspects of the glass-forming nature of the 
alloys: the thermodynamics, glass transition and struc- 
ture. There are three intrinsic limitations to these kinds 
of simulations: the interatomic potential, the system size 
and the timescale. We have reason to believe that the 
EMT interatomic potential is not a major limitation in 
this study. We have already discussed how much the 
physics of the binary amorphous alloy formation is based 
on size factors as well as the stability and structure of any 
nearby (in composition) intermetallic compounds. The 
fact that EMT parameters can be chosen to match quite 
closely the formation enthalpies of the two Mg-Cu inter- 
metallics means that the general bonding energetics are 
ious CNA-types vary smoothly with temperature. Fig.lHl reasonably well-represented. De Tendler et al.'*^ apphed 
shows how the positions and widths of these peaks vary. the empirical model of Miedema for alloy formation to 
One expects the widths to decrease as temperature de- compute the glass-forming region of the Mg-Cu system, 
creases, and this is indeed the case. Their heights in- The close agreement with experiment they found indi- 
crease, mostly to compensate for the narrowing: we have cates that there is nothing particularly unusual about 
already seen that the true measure of the weight of a this system. 

peak, the number of pairs associated with it, has only a This leads to the one feature of the Mg-Cu system 
small temperature dependence in the case of the second which is poorly described by our simulations: the ex- 
neighbor peaks. The splitting can now be seen as a natu- tent of the glass forming region. The width of the glass- 
ral consequence of this narrowing. It is also aided a little forming region is certainly a timescale-issue since the ac- 
by the decrease in weight of the 211 peak, which is in the cessible timescales preclude nucleation of a crystal phase 
middle, and the corresponding increase of the 333 peak more complex than fee; thus almost all compositions form 
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FIG. 12: Positions (upper panel) and widths (lower panel) 
of CNA components of second RDF peak for 50% Cu glass 
as a function of temperature. Squares: 333; diamonds: 211; 
triangles: 100. 

phases are locally very similar, the differences mostly be- 
ing those between "perfect" 555 pairs and "imperfect" 
544 and 433 pairs, and perfect and imperfect 333 and 
211 pairs respectively. 



E. Explanation of second peak splitting 

Our analysis indicates that the contributions from var- 
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a glassy phase upon cooling. Issues of length scale could 
conceivably also be relevant for the formation of the more 
complex Mg2Cu with its large unit cell. Of course, an 
advantage of being able to simulate glass-formation in a 
wide range of compositions is that it makes clearer that 
the splitting of the second peak and the glass transition 
are not coupled, since their dependences on temperature 
do not match. If one leaves aside crystallization, the fact 
that our results are largely independent of cooling rate, 
and the fact that the glass transition temperature for the 
eutectic composition matches the experimental one, sug- 
gest that the timescale is not otherwise a problem — until 
the onset of the glass transition itself of course; there 
the fast cooling rates lead to the broadening of the tran- 
sition compared to what would be expected experimen- 
tally. This of course does not rule out the possibility that 
there are relaxations that take place on time scales sig- 
nificantly longer than those of our slowest cooling rate, 
yet still fast enough to take place during the experimen- 
tal cooling. A possibility is that such relaxations might 
be associated with a length scale larger than our system 
size; thus our cooling rates are all slow enough to relax all 
structural rearrangements that are smaller than our sys- 
tem size, and thus we do not see any rate dependence, but 
perhaps we would see it in larger systems — the timescale 
issue and lengthscale issue would be in effect canceling 
each other out. However any such extra relaxations must 
be very low energy, because the residual enthalpies with 



respect to the crystal phases are as small as are mea- 
sured experimentally. Another "canceling" possibility is 
that defects in the interatomic potential, causing energy 
barriers to relaxation to be lower than they should, would 
lead to the relaxation times being lower and thus to the 
simulated cooling rate being more adequate than it other- 
wise should be. Guerdane and Teichlei»^ simulated Ni-Zr 
and ternary Ni-Zr-Al glass formation and obtained TgS 
higher than experimental ones by a few hundred Kelvin, 
which they explained as being due to the difference be- 
tween their cooling rate (lO^'^K/s) and the experimental 
one, which makes it surprising that we do not see such a 
discrepancy. If it is indeed due to too-low relaxation bar- 
riers, this may not matter so much for the purpose of ob- 
taining low temperature glassy configurations; however it 
may be relevant for future studies of plastic deformation. 
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